\chapter{Implementierung \label{implementierung}}
Dieses Kapitel präsentiert die Implementierungen diverser Filter inklusive der Funktionsplots der jeweiligen Kerne. Zunächst wird ein Testvolumen generiert, das als Referenz für alle Filter dient. Alle Implementierungen liegen als GLSL-Shader vor und werden durch ein \lstinline|#include| in den Shader \verb|mod_sampler3d.frag| von Voreen eingebunden. Dadurch können alle Rendering-Prozessoren, die diesen Shader einbinden, ebenfalls die neuen Filter nutzen. Dazu muss jeder Aufruf des gewrappten Hardwarefilters \lstinline|getVoxel()| durch den Funktionsaufruf des Filters ersetzt werden.
\section{Das Testvolumen \label{testvolumen}}
Um die Qualität der nachfolgend vorgestellten Filter einschätzen zu können, ist ein Vergleich zwischen den jeweiligen Renderings und einem Referenzrendering nötig. Das Referenzrendering des Testvolumens von S. R. Marschner und R. J. Lobb aus \cite[Abschnitt 6.1]{filter} wird mit dem Code aus Listing \ref{listingMarschnerLobb} berechnet. Es handelt sich um eine exakte Berechnung, bei der kein Filter zum Einsatz kommt. \\
Alle Filter rekonstruieren eine gesampelte Version des Testvolumens und werden in Kapitel \ref{benchmark} einem Vergleich mit diesem Referenzvolumen unterzogen.
\lstinputlisting[caption=Das Testvolumen von S. R. Marschner und R. J. Lobb,label=listingMarschnerLobb]{src/marschnerLobb.frag}
Das Ergebnis ist eine Wellenform (Abbildung \ref{marschnerLobb}), dessen hoher Kontrast - also die große Wertedifferenz zwischen benachbarten Samplepositionen - einen hohen Anspruch an das Renderingverfahren stellt und somit vorteilhaft bei der Beurteilung der Filterqualität ist.
\begin{figure}
\centering
\includegraphics[scale=\scale]{img/marschnerLobb.png}
\caption{Das Referenzvolumen von Marschner und Lobb}
\label{marschnerLobb}
\end{figure}
\section{Der Trilineare Filter \index{Trilinearer Filter}\label{trilinearFilter}}
Die simpelste Methode zur Rekonstruktion eines Zwischenwertes ist die lineare Interpolation zwischen den benachbarten Samples. Es ist das Standardverfahren zur Rekonstruktion, sodass es auch in der Hardware von Grafikkarten implementiert ist. Die dreidimensionale Variante der linearen Interpolation heißt \textit{trilineare} Interpolation. \\
In der Rekonstruktionstheorie entspricht die (tri-)lineare Interpolation einer Faltung mit der \textit{Dreiecksfunktion},  \index{Dreiecksfunktion} unter der Annahme, dass das Volumen in einem Einheitsgitter vorliegt. Da die Dreiecksfunktion (Abbildung \ref{dreiecksFunktion}) für Argumente mit einem Absolutwert größer eins verschwindet - respektive nur direkte Nachbarvoxel betrachtet werden - genügt es, für die diskrete Faltung den Radius auf eins zu beschränken.
\begin{figure}[h]
\centering
\includegraphics[scale=1]{img/dreiecksFunktion.pdf}
\caption{Die Dreiecksfunktion}
\label{dreiecksFunktion}
\end{figure}
\begin{figure}[h]
\includegraphics[scale=1]{img/trilinearSchema.pdf} 
\caption{Die trilineare Interpolation}
\label{figureTrilinear}
\end{figure} \\
Der lineare Filterkern erlaubt eine unkomplizierte Umsetzung, die auf eine Faltung und auf die in \ref{anwendung} beschriebene Separierung verzichtet: Wie in Abbildung \ref{figureTrilinear} zu sehen ist, berechnet sich die Interpolation eines Punktes (blau) zunächst durch lineare Interpolationen zwischen den umliegenden Voxeln (schwarz) entlang der x-Achse (a). Die resultierenden Interpolationen (rot) werden dann in Richtung der y-Achse ebenfalls linear interpoliert (b). Schließlich werden die neuen Interpolationen abermals linear in Richtung der z-Achse interpoliert (c), weswegen man insgesamt von einer trilinearen Interpolation spricht. \cite[Seite 23f]{real}\\ Listing \ref{listingTrilinear} zeigt eine Mögliche Implementierung des trilinearen Filters.
\lstinputlisting[caption=Die trilineare Interpolation,label=listingTrilinear]{src/trilinearFilter.frag}
Zunächst werden alle acht benachbarten Voxel bestimmt (Zeile 8 bis 16). Wie oben beschreiben wird zunächst in Zeile 24 bis 30 entlang der x-Achse, danach in Richtung der y-Achse (Zeile 34 bis 36), sowie in Zeile 39 an der z-Achse interpoliert. Das Ergebnis wird in den Rückgabewert geschrieben. \cite[Abschnitt 4.1]{filter} 
\section{Sinc Filter \index{Sinc Filter} \label{sincFilterFamilie}}
Wie in \ref{ideal} beschrieben, ist die $sinc$-Funktion der ideale Rekonstruktionsfilterkern. Da der $sinc$ andererseits relativ langsam gegen null konvergiert, erzwingt dies eine Faltung über einen dementsprechend großen Radius - was sich in der Praxis als problematisch erweist. Alternativ ist eine Modifikation des Sinc Filters möglich. Mittels einer Multiplikation der $sinc$ Funktion mit einer sogenannten \textit{Fensterfunktion}, \index{Fensterfunktion} durch die die Funktion schneller gegen null läuft, lässt sich die Faltung auf einen kleineren Radius beschränken. \\ Nachfolgend wird der Lanczos Filter und der abgeschnittene Sinc Filter vorgestellt.
\subsection{Der Lanczos Filter \index{Lanczos Filter} \label{lanczosFilter}} 
Der \textit{Lanczos Filter} ist ein möglicher gefensterter Sinc Filter. Als Fensterfunktion wird hierbei eine weitere $sinc$ Funktion genutzt:\[h(x) := lanczos(x) = sinc(\pi x) sinc(\frac{\pi x}{r})\] 
für einen festen Radius $r$. \cite{wikiLanczos}
\begin{figure}
\centering
\includegraphics[scale=1]{img/lanczosFunktion-3.pdf}
\caption{Die normierte Lanczos-Funktion für {$r = 3$}}
\label{lanczosFunktion}
\end{figure}
Abbildung \ref{lanczosFunktion} zeigt den normierten Lanczos für $r = 3$, der wie gewünscht schneller gegen null konvergiert als der Standard-Sinc Filter. \\
Eine simple Implementierung für $r = 2$ stellt sich wie folgt in Listing \ref{listingLanczosSimpel} dar:
\lstinputlisting[caption=Der Lanczos Filter (simpel) für {$r = 2$} ,label=listingLanczosSimpel]{src/simpleLanczosFilter.frag}
Ab Zeile 11 geschieht die eigentliche Faltung. Dazu wird zunächst in Zeile 15 die aktuelle Sampleposition von der normierten Form in eine Texturposition überführt. Zeile 17 bestimmt die dazu nächstgelegenen Voxel innerhalb des Radius. Die Differenz der Positionen bildet das Argument für den Lanczos Filter. Dieser berechnet sich ab Zeile 20 in einer weiteren Schleife, da für jede Dimension das Argument auf null überprüft werden muss. Schließlich wird das Signal \lstinline|getVoxel()|, multipliziert mit dem vorher berechneten Lanczos, aufsummiert und in den Alphakanal des Rückgabewerts geschrieben. (Zeile 27) Hierbei findet die in \ref{anwendung} beschriebene Separierung ihre Anwendung. 
\subsubsection{Optimierung}
Es fällt auf, dass der gleiche Lanczos mehrfach berechnet wird. Beispielsweise wird für die x-Achse der Lanczos bei jedem Schleifendurchlauf der beiden anderen Achsen jedes mal neu berechnet. Dieses Verhalten lässt sich durch die Vorausberechnung des Lanczos für jede Achse einzeln beheben.\\
Weniger offensichtlich ist, dass Zeile 22 des Listings \ref{listingLanczosSimpel} im Betrag konstant ist. Das Argument für den Sinus ist \[\pi (t - i) = \pi ( t - [n + t]) = \pi ( t - n - [t]) = \pi ( d - n)\] für $d := t - [t]$, $n := $ \lstinline|offset| und der Gaussklammer $[\cdot]$. Also unterscheidet sich das Argument jeweils um ein ganzzahliges Vielfaches von $\pi$. Wegen der Achsensymmetrie und der $\pi$-Periodizität des $sin$ gilt:\begin{eqnarray*}
sin(\pi(d-n)) &=& sin(\pi d - \pi n) \\
&=& - sin ( -\pi d + \pi n) \\
&=& (-1)^{n + 1} sin (- \pi d)\\
&=& (-1)^{n} sin (\pi d)
\end{eqnarray*}
Eine optimierte Version des Lanczos Filtes kann daher wie in Listing \ref{listingLanczos} aussehen.
\lstinputlisting[caption=Der Lanczos Filter (optimiert) für {$r=3$} ,label=listingLanczos]{src/lanczosFilter.frag}
Zeile 6 formt die Sampleposition wieder in Texturkoordinaten um und in Zeile 7 wird die nächstgelegene Texturkoordinate bestimmt. Dadurch lässt sich der Abstand \lstinline|distance| der beiden Positionen berechnen (Zeile 8), der fortan als Grundlage zur Auswertung des Lanczos dient. Ist die Sampleposition nahe genug an einem Voxel, wird die Berechnung abgebrochen und dieser Voxel als Ergebnis ausgegeben. Der gewählte Epsilonwert (Zeile 9) garantiert dabei, dass das Ergebnis zu über $ lanczos(0.001) = sinc(0.001 \pi) sinc(\frac{0.001 \pi}{3}) \approx 99.99 \% $ mit der Filterung des Lanczos übereinstimmt. Ab Zeile 16 werden schließlich alle benötigten Lanczoswerte für jede Dimension einzeln vorberechnet. Hierbei wird dem beschriebenen Umstand des für die speziellen Argumente alternierenden, aber im Betrag konstanten Sinus, durch eine angepasste Formel in Zeile 23 Rechnung getragen. Ab Zeile 32 wird die eigentliche Faltung mithilfe des vorberechneten Lanczos Filters durchgeführt. \\
Ein Vergleich der beiden Versionen des Lanczos zeigt, dass der nicht optimierte Filter $2(2r + 1)^3$ mal den Sinus berechnet und daher mit dem Radius $r$ als Eingabe in $\mathcal{O}(r^3)$ liegt, wohingegen der optimierte Lanczos nur $3(2r +1) +1$ mal den Sinus auswertet und damit nur lineare Kosten $\mathcal{O}(r)$ verursacht.\\
\\
Alle folgenden implementierten Filter verwenden dieselbe Struktur, wie die optimierte Version des Lanczos Filters.
\subsection{Der abgeschnittene Sinc Filter \index{Sinc Filter} \label{abgSincFilter}}
\begin{figure}
\centering
\includegraphics[scale=1]{img/sincFunktion.pdf}
\caption{Die normierte Sinc-Funktion}
\label{sincFunktion}
\end{figure}
Mit denselben Optimierungen wie denen des Lanczos Filters, reduzieren sich zwar alle Berechnungen des \textit{Sinc Filters} auf die Auswertung eines einzigen Sinus und auf Texture-Lookups, dennoch verhindert das kubische Wachstum der Anzahl von Lookups größere Radien im praktischen Einsatz. Daher ist eine Begrenzung auf kleine Radien wie zwei oder drei nötig. Die relativ großen Funktionswerte jenseits dieser Radien, die der Plot in Abbildung \ref{sincFunktion} offenlegt, lassen allerdings schon an dieser Stelle artefaktbehaftete Renderings vermuten. Das Ergebnis ist, wie alle anderen Renderings, im nächsten Kapitel abgebildet.
\lstinputlisting[caption=Der Sinc Filter für {$r = 3$} ,label=listingSinc]{src/sincFilter.frag}
Die Implementierung in Listing \ref{listingSinc} folgt dem bekannten Gerüst des Lanczos Filters (Listing \ref{listingLanczos}). In Zeile 17 wird der einzige erforderliche Sinus berechnet, der in Zeile 23 für die Vorberechnung des Filters benötigt wird. 
\section{Kubische Filter \index{Kubische Filter} \label{kubischeFilter}}
Die \textit{Kubischen Filter} bilden eine Schar von Splines. \index{Spline} In diesem Fall handelt es sich dabei um Polynome dritter Ordnung, die auf den Intervallen $(-\infty, -2], (-2,-1], (-1,1), [1,2)$ und $[2, \infty)$ definiert sind und zusätzlich an den Übergängen bis zur zweiten Ableitung übereinstimmen. \cite[Abschnitt 4.1]{filter} \cite[Seite 39]{num} Sie sind hier durch die Parameter $B$ und $C$ aus dem Intervall $[0, 1]$ gesetzt als\[h_{B,C}(x) := \frac{1}{6} \begin{cases} (12 - 9B - 6C)|x|^3 + \\ (-18 + 12B + 6C)|x|^2 + (6 - 2B) & \text{falls } |x| < 1 \\ (-B - 6C)|x|^3 + (6B + 30C)|x|^2 + \\ (-12B - 48C)|x| + (8B + 24C) & \text{falls } 1 \leq |x| < 2 \\
0 & \text{sonst}
 \end{cases}\]
\subsection{Der allgemeine Kubische Filter}
Die oben beschriebenen Splines lassen sich - wie in Listing \ref{listingkubisch} zu sehen - ihrer Definition entsprechend vollständig in Code umsetzen.
\lstinputlisting[caption=Der Kubische Filter ,label=listingkubisch]{src/cubicFilter.frag}
\subsection{Spezialfälle} 
Der B-Spline und der Catmull-Rom Filter stellen Spezialfälle des Kubischen Filters für fest gewählte Parameter $B$ und $C$ dar.

\pagebreak
\subsubsection{B-Spline Filter ($B = 1, C = 0$) \index{B-Spline Filter}
\label{bSplineFilter}} 
Für die Parameter $B = 1$ und $C = 0$ ergibt sich die in Abbildung \ref{bSplineFunktion} dargestellte zusammengesetzte Filterfunktion. Mit diesen Parametern bildet der Kubische Filter den sogenannten B-Spline. Die Implementierung des \textit{B-Spline Filters} (Listing \ref{listingBSpline}) ist dann ein Funktionsaufruf. 
\begin{figure}
\centering
\includegraphics[scale=1]{img/bSplineFunktion.pdf}
\caption{Die zusammengesetzte B-Spline Funktion}
\label{bSplineFunktion}
\end{figure}
\lstinputlisting[caption=Der B-Spline Filter ,label=listingBSpline]{src/bSplineFilter.frag}

\subsubsection{Catmull-Rom Filter ($B = 0, C = 0.5$) \index{Catmull-Rom Filter}
\label{catmullRomFilter}}
Abbildung \ref{catmullRomFunktion} zeigt den Filterkern des sogenannten \textit{Catmull-Rom Filters}. Er definiert sich durch die Parameterwahl $B = 0$ und $C = 0.5$ für die Kubischen Filter. Wie beim B-Spline Filter gestaltet sich die Implementierung (Listing \ref{listingCatmullRom}) sehr einfach.
\lstinputlisting[caption=Der Catmull-Rom Filter ,label=listingCatmullRom]{src/catmullRomFilter.frag}
\begin{figure}[]
\centering
\includegraphics[scale=1]{img/catmullRomFunktion.pdf}
\caption{Die zusammengesetzte Catmull-Rom Funktion}
\label{catmullRomFunktion}
\end{figure}
\section{Der Gauss Filter \index{Gauss Filter} \label{gaussFilter}}
Der \textit{Gauss Filter} ist aus der Bildbearbeitung bekannt. Dessen Filterkern ist dabei eine Glockenkurve beziehungsweise Normalverteilung. \[h(x) := gauss(x) = e^{-\frac{x^2}{2 \sigma^2}}\]
In Abbildung \ref{gaussFunktion} ist die Glockenkurve für $\sigma = 0.6$ zu sehen. Da dies eine typische Wahl für $\sigma$ ist und die weil die Funktion relativ schnell gegen null konvergiert, liefert der Algorithmus in Listing \ref{listingGauss}, der bereits ab einem Radius größer zwei abbricht, brauchbare Ergebnisse. \cite[Abschnitt 4.1]{filter}
\begin{figure}
\centering
\includegraphics[scale=1]{img/gaussFunktion-0,6.pdf}
\caption{Die Gauss-Funktion für $\sigma = 0.6$}
\label{gaussFunktion}
\end{figure}
\lstinputlisting[caption=Der Gauss Filter ,label=listingGauss]{src/gaussFilter.frag} 